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ABSTRACT 

We study the statistical properties of spherical harmonic modes of temperature 
maps of the cosmic microwave background. Unlike other studies, which focus 
mainly on properties of the amplitudes of these modes, we look instead at their 
phases. In particular, we present a simple measure of phase correlation that 
can be diagnostic of departures from the standard assumption that primor- 
dial density fluctuations constitute a statistically homogeneous and isotropic 
Gaussian random field, which should possess phases that are uniformly ran- 
dom on the unit circle. The method we discuss checks for the uniformity of 
the distribution of phase angles using a non-parametric descriptor based on 
the use order statistics, which is known as Kuiper's statistic. The particular 
advantage of the method we present is that, when coupled to the judicious use 
of Monte Carlo simulations, it can deliver very interesting results from small 
data samples. In particular, it is useful for studying the properties of spherical 
harmonics at low I for which there are only small number of independent values 
of m and which therefore furnish only a small number of phases for analysis. 
We apply the method to the COBE-DMR and WMAP sky maps, and find de- 
partures from uniformity in both. In the case of WMAP, our results probably 
reflect Galactic contamination or the known variation of signal-to-noise across 
the sky rather than primordial non-Gaussianity. 

Key words: cosmic microwave background - cosmology: theory - large-scale 
structure of the Universe - methods: statistical 



1 INTRODUCTION 

A crucial ingredient of many cosmological models is the idea that galaxies and large-scale structure in the Universe 
grew by a process of gravitational instability from small initial perturbations. In the most successful versions of this 
basic idea, the primordial fluctuations that seeded this process were generated during a period of inflation which, 
in its simplest form, is expected to produce fluctuations with relatively simple statistical properties (Starobinsky 
1979, 1980, 1982; Guth 1980; Guth & Pi 1981; Linde 1982; Albrecht & Steinhardt 1982). In particular the primordial 
density field in these models is taken to form a statistically homogeneous (i.e. stationary) Gaussian random field 
(Bardeen et al. 1986). This basic paradigm for structure formation has survived numerous observational challenges, 
and has emerged even stronger after recent confrontations with the 2dF Galaxy Redshift Survey (2dFGRS; Percival 
et al. 2001) and the Wilkinson Microwave Anisotropy Probe (WMAP; Hinshaw et al. 2003). So successful has the 
standard paradigm now become that many regard the future of cosmology as being largely concerned with improving 
estimates of the parameters of this basic model rather than searching for alternatives. 

There are, however, a number of suggestions that this confidence in the model may be misplaced and the focus 
on parameter estimation may be somewhat premature. For example, the WMAP data have a number of unusual 
properties that are not yet completely understood (Efstathiou 2003; Chiang et al. 2003; Dineen & Coles 2003; 
Eriksen et al. 2003). Among the possibilities suggested by these anomalies is that the cosmic microwave background 
(CMB) sky is not statistically homogeneous and isotropic, and perhaps not Gaussian either. The latter possibility 
would be particularly interesting as it might provide indications of departures from the simplest versions of inflation 
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(e.g. Linde & Mukhanov 1997; Contaldi, Bean & Magueijo 2000; Martin, Riazuelo & Sakellariadou 2000; Gangui, 
Pogosian & Winitzki 2001; Gupta et al. 2002; Gangui, Martin & Sakellariadou 2002; Bartolo, Matarrese & Riotto 
2002). Whether WMAP is hinting at something primordial or whether there are systematic problems with the data or 
its interpretation is unclear. Either way, these suggestions, as well as general considerations of the nature of scientific 
method, suggest that it is a good time to stand back from the prevailing paradigm and look for new methods of 
analysis that are capable of testing the core assumptions behind it rather than taking them for granted. 

In this paper we introduce a new method for the statistical analysis of all-sky CMB maps which is complementary 
to usual approaches based on the power-spectrum but which also furnishes a simple and direct test of the statistical 
assumptions underpinning the standard cosmological models. Our method is based on the properties of the phases 
of the (complex) coefficients obtained from a spherical harmonic expansion of all-sky maps. The advantages of this 
approach are that it hits at the heart of the "random-phase" assumption essential to the definition of statistically 
homogeneous and isotropic Gaussian random fields. Perhaps most importantly from a methodological point of view, 
is intrinsically non-parametric and consequently makes minimal assumptions about the data. 

The layout of this paper is as follows. In the next Section we discuss some technical issues relating to the properties 
of spherical harmonic phases which are necessary for an understanding of the analysis we present. In Section 3 we 
explain a practical procedure for assessing the presence of a particular form of departure from the random phase 
hypothesis using Kuiper's statistic. In Section 4 we discuss results obtained by applying the method to COBE-DMR 
maps and data from WMAP, as well as some toy examples. We discuss the outcomes and outline ideas for future 
work in Section 5. 



2 TECHNICAL PREAMBLE 
2.1 Fourier Phases 

The method we discuss in this paper is based on recent results arising from the study of correlations between Fourier 
domain for random fields 5(x) defined over flat two- or three-dimensional spaces. In two dimensions, which is the 
case closest to the spherical expansion we use in this paper, the Fourier expansion is of the form 

6{x) = 6{x,y) = ^^6{k)exp{ik-x), (1) 
with k = (k x ,k y ) and the Fourier transform <5(k) is complex, i.e. 

5(k) = |5(k)|exp[#(fc x ,fe y )]. (2) 

The quantity 0(k) is the phase of the mode corresponding to wavevector k. A two-dimensional flat surface is a 
reasonable approximation to a small patch of the sky, so a Fourier approach has been taken in studies of high- 
resolution CMB maps (Bond & Efstathiou 1987; Coles & Barrow 1987). The extension to three-dimensions is trivial. 

In orthodox cosmologies, initial density fluctuations constitute a statistically homogeneous and isotropic Gaussian 
random field (Bardeen et al. 1986). Strictly speaking this means that the real and imaginary parts of <5(k) are drawn 
independently from Gaussian distributions with a variance that depends only on |k|. If this assumption holds, the 
phases 0(k) are independent and uniformly random on the interval [0, 2tt]; for more details see Watts & Coles (2003). 
On the other hand, the central limit theorem comes into play whenever a large number of independent Fourier modes 
are involved in the superposition represented by Equation (1). As a result, the random-phase hypothesis on its own 
suffices to define a weak form of Gaussianity (Bardeen et al. 1986). As long as the evolution of density perturbations 
is linear, initially random phases will remain so. But in the non-linear regime, mode-mode interactions can introduce 
non-randomness into the distribution of phases. Characterizing non-linear clustering growth using phase information 
has been a challenge for a number of years (Peebles 1980; Ryden & Gramman 1991; Scherrer, Melott & Shandarin 
1991; Soda & Suto 1992; Jain & Bertschinger 1996, 1998), but recently there have been some important advances 
in developing practical methods based on the use of Fourier phases (Chiang & Coles 2000; Coles & Chiang 2000; 
Chiang 2001; Chiang, Coles & Naselsky 2002; Watts, Coles & Melott 2003; Matsubara 2003; Hikage, Matsubara & 
Suto 2003) and also on relating phase information to more traditional statistical approaches based on quantities such 
as the bispectrum (Watts & Coles 2003) . 

One of the difficulties encountered when constructing methods based on phase information is how to construct 
quantitative measures of association appropriate for a circular measure such as 0(k). Chiang & Coles (2000) pointed 
out that traditional measures of covariance of the form (XY) between variates X and Y do not work for two phases 
0(ki) and (f)(k.2) because of the circular nature of phase angles. They also considered the obvious generalization to 
expectations of the form 
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(exp[i0(ki)] exp[i0(k 2 )]} = (exp[i(fa + fa)]), (3) 

for any two wavevectors ki and k 2 , except for the case ki = — k 2 . Such correlations are zero for any random field 
that is statistically homogeneous and isotropic, i.e. one that has statistical properties that are both translation and 
rotation invariant, whether or not it is Gaussian. This essentially means that the sum of phases fa + fa is also 
random; a similar argument can be made for phase differences fa —fa. The case ki = — k 2 is excluded from the 
above argument because the requirement that the original field 8 is real imposes the constraint that 

5(k) = S(-k)* (4) 

so that 0(k) = — fa— k), fa + fa is identically zero and the expectation value in Equation (3) is unity. 

One of the useful tricks developed by Chiang & Coles (2000) for the analysis of numerical iV-body simulations was 
to consider phase differences instead of the phases themselves. Such simulations typically employ periodic boundaries 
and a finite grid, so the sum in equation (5) is taken over equally-spaced discrete wavenumbers up to some maximum 
value determined by the size of the periodic box. Taking the spacing to be unity one can most simply construct phase 
differences along one direction in k-space, such as 

D k = fak + 1) - fak) = arg[<5(fc + 1)5* (k)]; (5) 

a straightforward generalization of this approach leads to a technique known as phase mapping (Chiang, Coles & 
Naselsky 2002; Chiang, Naselsky & Coles 2002) but we shall stick with the simple form given in equation (5). On top 
of its simplicity this approach has a number of advantages. One is that although translating the coordinate system by 
an amount x shifts the phases by kx, the phase difference changes by the same amount for each wavevector. Secondly, 
if the phases are uniformly random on the interval [0, 2ir], then the phase differences are also uniformly random on 
the same interval. Finally, the behaviour of the phase differences is relatively simple to visualize, and interpret in 
terms of the dynamics of non-linearly evolving structure (Coles & Chiang 2000). 

However the preceding argument seems to demonstrate that statistical homogeneity and isotropy require that 
< expiDfe >= or, in other words, that the phase differences are random. So how can D k be useful statistically? The 
answer to this question is that although the expectation value taken over the ensemble of possible configurations of a 
random field may indeed be zero, any average taken over a finite region within single realization will not correspond 
to this value. Departures from the global expectation value contain information about higher-order fluctuations. In 
the case of non-linear gravitational clustering developing within a finite numerical box, what is key is the scale of 
the non-linearity compared with the size of the box. If non-linearity exists only on small scales, then the D k are very 
close to uniform on the unit circle and the box average is close to the ensemble average ("the box is a fair sample"). 
If the scale of structure is close to the scale of the box then D k becomes extremely non-uniform (Coles & Chiang 
2000). Using a bigger box actually makes it harder to characterize non- linear structure on a given scale using phase 
differences. For related comments see Hikage, Matsubara & Suto (2003). 



2.2 Spherical Harmonic Phases 

We can describe the distribution of fluctuations in the microwave background over the celestial sphere using a sum 
over a set of spherical harmonics: 

— oo m — 

A(6, fa = ^0A)-T) ^^(0, (6) 

£ — 1 m — — I 

Here A(9, fa is the departure of the temperature from the average at angular position (9, fa on the celestial sphere in 
some coordinate system, usually Galactic. The Yi m (6,fa are spherical harmonic functions which we define in terms 
of the Legendre polynomials P im using 



Y lm (6,fa = (-ir^ {2l +^ m ™ )[ P lm (cos8)e X p(irnfa, (7) 

i.e. we use the Condon- Shortley phase convention. In Equation (6), the a !im are complex coefficients which can be 
written 

ai, m = \ai, m \exp[i<i>i, m ]. (8) 

Note that, since A is real, the definitions (7) and (8) require the following relations between the real and imaginary 
parts of the aj jm : if m is odd then 

SR(oi, m ) = -5R(oi _ m ), S(a iim ) = 9(a; _ m ); (9) 
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while if m is even 



»(oj, m ) = »(o,, 



— m 



9(a;, m ) 



3(ai,-m); 



(10) 



and if m is zero then 



9?(oi,m) = 0. 



(11) 



From this it is clear that the m — mode always has zero phase, and there are consequently only I independent phase 
angles describing the harmonic modes at a given I. 

If the primordial density fluctuations form a Gaussian random field as described in Sec. 2.1, then the temperature 
variations induced across the sky form a Gaussian random field over the celestial sphere. By analogy with the above 
discussion above, this means that 

(oi, m Oi* jm /) = Ci5u>6 mrn ' , (12) 

where Ci is the angular power spectrum, the subject of much scrutiny in the context of the cosmic microwave 
background (e.g. Hinshaw et al. 2003), and 5 xx i is the Kronecker delta function. Since the phases are random, the 
stochastic properties of a statistically homogeneous and isotropic Gaussian random field are fully specified by the Ci, 
which determines the variance of the real and imaginary parts of a;, m both of which are Gaussian. 

To look for signs of phase correlations in CMB maps one can begin with a straightforward generalization of the 
ideas introduced above for Fourier phases. There are, however, some technical points to be mentioned. First, the set of 
spherical harmonics, and consequently the phases of the harmonic coefficients, is defined with respect to a particular 
coordinate system. This is also true in the Fourier domain but, as we discussed above, dealing with translations of 
the coordinate system in that case is relatively straightforward because k x and k y are equivalent. This is not so for 
spherical harmonics: the definitions of I and m are quite different, and the distinction is introduced when the «-axis 
of a three-dimensional coordinate system is chosen to be the polar axis. Changing this axis muddles up the spherical 
harmonic coefficients in a much more complex way than is the case for Fourier modes. 

Another point worth stressing is that the sphere is a finite space, so when analyzing the statistical properties of 
maps on a single sky one never has an exact ensemble average. This is less of an issue when dealing with spherical 
harmonic modes at high /, because the average over a single sky is then close to that over the probability distribution, 
but at low I there is the perennial problem of so-called "cosmic variance". Suppose one were to construct phase 
differences in the manner of Section 2.1, i.e. 

D m (l) = (f>l, m +l — (f>l, m - (13) 

Then the distribution of the D m for high I should tend to uniform if the field is statistically homogeneous over the 
sphere whether or not the temperature fluctuations are Gaussian. However, at low I, the finiteness of the sky come 
into play and the distribution over the small number of modes available will depart from the uniform case. These 
fluctuations at low I depend on higher-order statistical information about the form of the fluctuations, so in this case 
cosmic variance can actually be helpful. 

We shall explain how we cope with this aspect of phases and the behaviour in more detail in Section 3. For the 
time being we note that, if the orthodox cosmological interpretation of temperature fluctuations is correct, the phases 
of the a; iJrl should be random and so should phase differences of the form <f>i >m+ i — 4>i >rn and 0; + i jm — (j>i, m - The aim 
of this paper is introduce a method of checking whether this is so. 

2.3 Departures from Standard Statistics 

In studying the properties of phase correlations of CMB maps, our intention is to characterize departures from the 
standard cosmological framework. The usual hypothesis, that primordial density fluctuations form a statistically ho- 
mogeneous and isotropic Gaussian random field, results in uncorrelated phases. Most discussions of phase correlations 
in the literature refer to them as diagnostic of non-Gaussianity in some form. What this actually means is usually 
that the field in question is not a homogeneous and isotropic Gaussian random field. Fields can be defined which are 
non-Gaussian but are both homogeneous and isotropic; Coles & Barrow (1987) give examples. Such fields certainly 
do not have random phases, but the form of phase association that characterizes them can be quite complex (Watts 
& Coles 2003; Matsubara 2003; Hikage, Matsubara & Suto 2003). However, Watts & Coles (2003) established an 
important connection between the form of non-Gaussianity exhibited by such fields and the bispectrum. The bispec- 
trum is now a standard part of the cosmologists armoury for studies of large-scale structure in three dimensions (e.g. 
Peebles 1980; Luo 1994; Heavens 1998; Matarrese, Verde & Heavens 1997; Scoccimarro et al. 1998; Scoccimarro et 
al. 1998; Scoccimmarro, Couchman & Frieman 1999; Verde et al. 2000, 2001, 2002) and for CMB studies (Ferreira, 
Magueijo & Gorski 1998; Heavens 1998; Magueijo 2000; Contaldi et al. 2000; Phillips & Kogut 2001; Sandvik & 
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Magueijo 2001; Santos et al. 2001; Komatsu et al. 2002; Komatsu et al. 2003). The bispectrum, being essentially 
the three-point correlation function in harmonic space, can be generalized to higher order polyspectra (Stirling & 
Peacock 1996; Verde & Heavens 2001; Hu 2001; Cooray 2001). However the estimation of these polyspectra generally 
involves taking averages that assume that the field one is dealing with is indeed statistically homogeneous and they 
are therefore not in themselves tests of that hypothesis. 

On the other hand, it is also possible for fields to be Gaussian in some sense, but not strictly Gaussian in 
the sense we defined in Sec. 2.1. As a consequence of strict Gaussianity, all the n-dimensional joint probability 
distributions of the field values at n different spatial locations are multivariate Gaussian (Bardeen et al. 1986). Fields 
can be constructed in which, say, the one-point distribution is Gaussian but the higher-order joint distributions 
are not. These would be Gaussian in some sense, but again would not have random phases and may or may not 
be statistically homogeneous and isotropic. In the Fourier or spherical harmonic domain, one can construct fields in 
which the harmonic coefficients are Gaussian but not independent. Such fields are Gaussian, but either statistically 
inhomogeneous or statistically anisotropic and would have non-random phases. 

Any or all of these deviations from simple Gaussianity might be expected to occur in CMB in various circum- 
stances. Statistical inhomogeneity and/or anisotropy over the sky might result from incorrectly subtracted foreground 
emission from the Galaxy. Alternatively, the scanning pattern used to produce all-sky map will generally not sample 
the sky uniformly resulting in a signal-to-noise that varies over the celestial sphere. Perhaps more exotically, compact 
universe models with a non-trivial topology should result in repeated features in the temperature pattern on the sky 
(Levin, Scannapieco & Silk 1998; Scannapieco, Levin & Silk 1999; Levin 2002; Rocha et al. 2002). All these possibili- 
ties should result in some form of phase correlation, and indeed there is already some evidence from the preliminary 
release of the WMAP data (Chiang et al. 2003; Eriksen et al. 2003). 

The focus of most discussions of higher-order CMB statistics is the possible detection of primordial non- 
Gaussianity. As we discussed in the introduction, most versions of the inflationary universe idea produce primordial 
fluctuations that are extremely close to the strictly Gaussian form. However, as we mentioned in the Introduction, 
there are circumstances in which the primordial fluctuations could be non-Gaussian. In such cases there will be 
information concerning the level of non-Gaussianity in the distribution of phases, but this usually appears in a form 
that is picked up by the bispectrum (Watts & Coles 2003). 

Clearly there are many different possible forms of departure from the standard statistical assumption and conse- 
quently many possible signatures in the distribution of phases. Our aim in this paper is to look for a simple method 
that is a powerful complement to the usual approach via the polyspectra and hopefully sheds some light on previous 
debates about departures from Gaussianity in the COBE-DMR data (Ferreira, Magueijo & Gorski 1998; Pando, Valls- 
Gabaud & Fang 1998; Bromley & Tegmark 1999; Magueijo 2000). As we shall see, a test based on the distribution of 
harmonic phases and/or phase differences fits the bill rather nicely. 



3 TESTING FOR PHASE CORRELATIONS 
3.1 Kuiper's Statistic 

The approach we take in this paper is to assume that we have available a set of phases <^ im corresponding to a set 
of spherical harmonic coefficients a;. m obtained from a data set, either real or simulated. We can also form from 
these phases a set of phase differences as described in the previous section. Let us assume, therefore, that we have n 
generic angles, 0i . . . n . Under the standard statistical assumption these should be random, apart from the constraints 
described in Section 2.2. The first thing we need is a way of testing whether a given set of phase angles is consistent 
with being drawn from uniform distribution on the unit circle. This is not quite as simple as it seems, particularly 
if one does not want to assume any particular form for actual distribution of angles, such as a bias in a particular 
direction; see Fisher (1993). Fortunately, however, there is a fully non-parametric method available, based on the 
theory of order statistics, and known as as Kuiper's statistic (Kuiper 1960). 

Kuiper's method revolves around the construction of a statistic, V, obtained from the data via the following 
prescription. First the angles are sorted into ascending order, to give the set {9i, . . . , 9 n }. It does not matter whether 
the angles are defined to lie in [0, 2n], [— 7T, +ir] or whatever. Each angle 9i is divided by 2-7T to give a set of variables 
Xi, where i = 1 . . . n. From the set of Xi we derive two values and S~ where 




(14) 



and 




(15) 
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Kuiper's statistic, V, is then denned as 

V = {S++S~)- ( ^ + 0.155 + ^ ) . (16) 



Anomalously large values of V indicate a distribution that is more clumped than a uniformly random distribution, 
while low values mean that angles are more regular. The test statistic is normalized by the number of variates, n, in 
such a way that standard tables can be constructed to determine significance levels for any departure from uniformity; 
see Fisher (1993). In this context, however, it is more convenient to determine significance levels using Monte Carlo 
simulations of the "null" hypothesis of random phases. This is partly because of the large number of samples available 
for test, but also because we can use them to make the test more general. 

The first point to mention is that a given set of phases, say belonging to the modes at fixed I is not strictly speaking 
random anyway, because of the constraints noted in equations (9)-(ll). One could deal with this by discarding the 
conjugate phases, thus reducing the number of data points, but there is no need to do this when one can instead 
build the required symmetries into the Monte Carlo generator. 

In addition, suppose the phases of the temperature field over the celestial sphere were indeed random, but 
observations were available only over apart of the sky, such as when a Galactic cut is applied to remove parts of 
the map contaminated by foregrounds. In this case the mask may introduce phase correlations into the observations 
so the correct null hypothesis would be more complicated than simple uniform randomness. As long as any such 
selection effect were known, it could be built into the Monte Carlo simulation. One would then need to determine 
whether V from an observed sky is consistent with having been drawn from the set of values of V generated over the 
Monte Carlo ensemble. 

There is also a more fundamental problem in applying this test to spherical harmonic phases. This is that a 
given set of a Jjm depends on the choice of a particular coordinate axis. A given sky could actually generate an infinite 
number of different sets of <j)i,m because the phase angles are not rotationally invariant. One has to be sure to take 
different choices of «-axis into consideration when assessing significance levels, as a random phase distribution has 
no preferred axis while systematic artifacts may. A positive detection of non-randomness may result from a chance 
alignment of features with a particular coordinate axis in the real sky unless this is factored into the Monte Carlo 
simulations to. For both the real sky and the Monte Carlo skies we therefore need not a single value of V but a 
distribution of F-values obtained by rotating the sky over all possible angles. A similar approach is taken by Hansen, 
Marinucci & Vittorio (2003). This method may seem somewhat clumsy, but a test is to be sensitive to departures 
from statistical homogeneity one should not base the test on measures that are rotationally invariant, such as those 
suggested by Ferreira, Mageuijo & Gorski (1998) as these involve averaging over the very fluctuations one is trying 
to detect. 



3.2 Rotating the a;, m 

In view of the preceding discussion we need to know how to transform a given set of a;, m into a new set when the 
coordinate system is rotated into a different orientation. The method is fairly standard, but we outline it here to 
facilitate implementation of our approach. 

Any rotation of the cartesian coordinate system S{x,y,z} i— » S'{x,y, z} can be described using a set of three 
Euler angles a, f3, 7, which define the magnitude of successive rotations about the coordinate axes. In terms of a 
rotation operator D(a,/3, 7), defined so that a field f(r,9,(p) transforms according to 

D(a, /3, 7 )/(r, 6, 4>) = f (r, 0, 4>) = f{r, 9', 0') , (17) 
a vector r is transformed as 

r' = D(0, 0, 7)D(0, f3, 0)D(a, 0, 0)r = D(a, /3, 7)1- . (18) 
Here D is a matrix representing the operator D, i.e. 



D(a,/3, 7) = —sin 7 cos 7 1 —sin a cos a . (19) 




The Wigner D functions describe the rotation operator used to realise the transformations of covariant components 
of tensors with arbitrary rank of rank I. The functions, written as D m ,, transform a tensor from S{x,y,z} to 
S'{x' , y' , z'}. Consider a tensor F( >m (#, <f>) defined under the coordinate system S and apply the rotation operator we 
get: 
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Figure 1. Mapping the symmetries (see text). 



D(a, 0, y)Y,, m (0, 0) = Y Um , (#', <f>') = ^ Y,, m (0, 4>)D l m , m , (a, (3, 7 ) 



(20) 



This means that the transformation of the tensor under the rotation of the coordinate system can be represented as 
a matrix multiplication. An alternative representation of the D functions D l m m , in terms of the Euler angles is 



D l m , m ,(a,(3, 1 )=e- lma d l m , m ,((3)e- 



(21) 



In this representation the parts of the D function depending upon the individual Euler angles are separated. The 
d l m m , form the elements of a rotation matrix and only depend upon the angle (3. Using this representation makes 
calculating the values of the D m m , easier since the d l mm i possess a number of symmetry properties meaning that 
although one quarter of the elements need to be directly calculated the values of the remainder can be inferred. Figure 
(1) shows how the following symmetries can be used to map the elements in the top quadrant of the rotation matrix 
into each of the other quadrants. First, 



(-1) d m 



which maps elements in region I into region II as shown in Figure 1. Next, 
which maps elements in region I into region III. And finally, 

il / -i\m— m -jl 

^*"tn,m' V *-) U' — m. — m'' 



(22) 
(23) 
(24) 



This maps elements in region I into region IV. The relations (22) to (24) can be used in conjunction with Equation 
(21), to find all the elements of the Wigner D functions from the values of d l m m , in the upper quadrant. Similar, but 
more complicated, symmetry relations exist for the D l m m , . 

For a given value of I we can construct a sum of spherical harmonics in the coordinate system S'{r, 6' , 4>'} as 
follows 



Yl =^2a'i,m'Yi,m'(0' A')- 



(25) 



We can use this and the transformation in equation (20) to find the transformed sum, Yj in the coordinate system 
S{r,9^} 



m' \ m / 



Q>l,m Yl,m{9, 4>) 



m \ m 



This gives an expression for the rotated coefficients 
]ai,m'D l mtm ,(a,/3,j). 



(26) 



(27) 



Finding the rotated coefficients therefore requires a simple matrix multiplication once the appropriate D function is 
known. To apply this in practice one needs a fast and accurate way of generating the matrix elements D l m m , for the 
rotation matrix. There are (21 + l) 2 elements needed to describe the rotation of each mode and the value of each 
element depends upon the particular values of (a, /3, 7)used for that rotation. 
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Figure 2. The basic steps involved in the iteration method described in the text for increasing matrix size: {left) 1=1, (middle) 
1 = 2, and (right) 1=3. 



Varshalonich, Moskalev & Khersonskii (1988) provide a set of matrices for the d l mrn ,((3) for the modes / 
1/2, 1, 3/2, . . . , 5. For example, the case I — 1 is 




(28) 



Fortunately, rather than working these matrices out by hand, Varshalovich et al. (1988) includes a set of iteration 
formulae, which, in conjunction with the symmetry relations for the d l m m ,, can be used to generate matrices for the 
higher modes from the previous ones. The first iteration formula provides a way of generating an element in the mode 
I + 1 by using the equivalent element in the matrices for the modes I and I — 1. 



D 



= f{D l m , 
= (Z + 1)(2Z + 1 



V(l 2 -mm 2 -m^)D^ m , + (cos/3 - £>;„. , 



v /((/ + l)2_ m 2)(( Z+1 )2_ m ,2) 



(29) 



Applying this formula repeatedly allows the central nine elements to be generated for all of the matrices. 

Another iteration formula uses the values already in the matrix to calculate the element directly to the left. For 
example the element Dq 2 could be generated using the elements and Dq . 



f ( n l , n l , \ 

X 

y/(l-m')(l + m' + 1) 
(— m + m' cos f3)D l m r 
sin (3 



m' + !)£>; 



(30) 



Similar iteration functions exist to move up a matrix and across a matrix to the right. 

Figure (2) shows the iteration scheme. Iteration formula (29) is used to generate the elements coloured in yellow 
from the elements coloured in cyan. The elements coloured in yellow are then used to calculate the elements coloured 
in blue. The yellow and blue elements are used to calculate the elements coloured in red. The red elements are used 
to calculate the elements coloured in green. The symmetry relations are used to calculate the elements in the other 
three quadrants. Once the I = 3 matrix has been calculated it can be used with the I = 2 matrix to find the fourth 
matrix, and so on and so forth. 

This method is fairly straightforward to implement but care is needed. The main problem with it is that the 
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second iteration formula shown, equation (30) includes a term in which a factor is divided by sin (3. If (3 ~ or (3 ~ tv 
then sin (3 ~ leading to numerical problems. 

An alternative approach is to calculate each d l mm , using factorials. The procedure encoded the definition of the 
d m m , given in Edmonds (1960): 



d m , m > (0) = 



{I + m'y.jl - m')\ 
(l + m)\(l-m)\ 



E 



! + m 
rri — k 



I — rn 
k 



(-1)' 



x (cos/3) 



2fc + m +m 



(sin (3) 



ll-lk-m 1 ' - 



(31) 



In this expression the variable k varies between and the smallest of I — m and I — m' . This method avoids the 
problem of dividing through by (sin/3) n when (3 is small. The results are far more stable and satisfy both of the 
techniques given above for testing the rotation code. The procedure can be tested using a random number generator 
to produce matrices for the full range of possible Euler angles. 

There are two main tests that can be made of the results. First, the rotation code should generate the correct 
values for the d l m m , matrices given in Varshalovich et al. (1988). Since the iteration formula are for the D l m m , they 
will be applicable to the d l m m , if a and 7 are set to zero. Choosing a value of (3 and using it in matrices (25) and (26) 
should reproduce the precise calculation of the I = 3, 4 and 5 matrices as given by Varshalovich et al. (1988). Second, 
the calculated matrices should be rotation matrices, so they should preserve the length of a vector. We define the 
vector of coefficients Ai for a particular mode I as 



Ai — (di,-i, a^-i+i, . . . , Oi,_i, a;,o, a;,!, . . . a(,i-i, a;,;) 
and use this to define the length of A t to be 

\A,\ 




(32) 



(33) 



The length should be preserved under the transformation. This allows the rotation code to be tested for any values 
of I > 5 and for 0,7/0. 



3.3 Implementation 

In order to apply these ideas to make a test of CMB fluctuations, we first need a temperature map from which we can 
obtain a measured set of ai. m . Employing equation 27 with some choice of Euler angles yields a rotated set of the Oj, m . 
It is straightforward to choose a set of angles such that random orientations of the coordinate axis can be generated. 
The Wigner D function needed to compute this transformation is generated using the techniques described in Section 
3.3. For the purposes of this paper we only compute the effect of rotation on low values of I. There is no difficulty 
in principle in extending this method to very high spherical harmonics, but the computational generation of Wigner 
D scales by ~ I 2 so have chosen to limit ourselves to £=20. Once a rotated set has been obtained, Kuiper's statistic 
is calculated from the relevant transformed set of phases. For each CMB map, 3000 rotated sets are calculated by 
this kind of resampling of the original data, producing 3000 values of V cm b. The values of the statistic are binned to 
form a measured (re-sampled) distribution of V cm b- The same procedure is applied to the 1000 Monte Carlo sets of 
a; jm drawn from a uniformly random distribution, i.e. each set is rotated 3000 times and a distribution of Vmc under 
the null hypothesis is produced. These realizations are then binned to created an overall global average distribution 
under the null hypothesis. In the following applications, the distribution of the 3000 values of V was split into 100 
equally-spaced bins ranging from V = to 3.0 for different I and 100 equally-spaced bins ranging from V = to 5.2 
for different m. 

In order to determine whether the distribution of V cmb is compatible with a distribution drawn from a sky with 
random phases, we use a simple % 2 test, using 

X 2 = £ (34) 

where the summation is over all the bins and /; is the number expected in the ith bin from the overall average 
distribution. The larger the value of x 2 the less likely the distribution functions are to be drawn from the same parent 
distribution. Values of Xmc are calculated for the 1000 Monte Carlo distributions and Xcmb ls calculated from the 
distribution of Kmb- If the value of Xcmb i s greater than a fraction p of the values of Xmc> then the phases depart 
from a uniform distribution at significance level p. We have chosen 95 per cent as an appropriate level for the level 
at which the data are said to display signatures that are not characteristic of a statistically homogeneous Gaussian 
random field. 
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Table 1. Results for non-Gaussian skies based on phases between consecutive values of m at a fixed I. The results show the 
significance levels of detected departures from random phases for simulated distributions corresponding to (i) coupled phases, 
(ii) a Cardioid distribution and (iii) a wrapped Cauchy distribution. In the latter two cases various parameter choices are shown 
as described in the text. These results arc obtained from 1000 Monte Carlo skies. 



4 APPLICATIONS 
4.1 Toy Models 

As a first check on the suitability of our method for detecting evidence of non-Gaussianity in CMB skymaps, sets 
of spherical harmonic coefficients were created with phases drawn from non-uniform distributions. Three types of 
distributions were used to test the method: (i) the phases were coupled (ii) the phases were drawn from a cardioid 
distribution and (iii) the phases were drawn from a wrapped Cauchy distribution; see Fisher (1993). The fake sets 
of coefficients were created in a manner that insured that the orthonormality of the coefficients were preserved as 
outlined in equations (9)-(ll). 

The first non-uniform distribution of phases is constructed using the following relationship 

0m = 4>m-l + k (35) 

where <$> m= -i is a random number chosen between and 2n. The results of the simulation with k = 1 are shown 
in panel (i) of Tables 1 and 2. The sky is taken to be non-Gaussian if Xsky ls larger than 95 % of Xmc f° r a 
particular mode. The phase difference method should be particularly suited to detecting this kind of deviation from 
Gaussianity because the phase differences are all equal in this case; the distribution of phase differences is therefore 
highly concentrated rather than uniform. Nevertheless, the non- uniformity is so clear that it is evident from the 
phases themselves, even for low modes. On the other hand, assessing non- uniformity from the quadrupole mode I — 2 
on its own is difficult even in this case because of the small number of independent variables available. 

The cardioid distribution was chosen because studies of phase evolution in the non-linear regime have shown 
phases rapidly move away from their initial values, however, they wrap around many multiples of 2-7T and the 'observed' 
phases appear random (Chiang and Coles 2000; Watts et al. 2003). This evolution can be distinguished by the 
phases differences which appear to be drawn from a roughly cardioid distribution (Watts et al. 2003). The cardioid 
distribution has a probability density function given by 

f(0) = -!-[l + 2 / 9cos(0)], < 9 < 2tt,0 < p < 1/2, (36) 

where p is the mean resultant length. As p — > the distribution converges to a uniform distribution. The kurtosis of 
the distribution is given by 

(1-P) 2 

Sets of coefficients were produced from distributions with p = 0, 0.2, 0.4, 0.5. The distributions themselves were 
generated using rejection methods. The results are given in panel (ii) of Tables 1 and 2. The results with p = are 
not shown: they indicate the coefficients are taken from a Gaussian random field, as expected. Furthermore, with 
p — 0.5 the non-uniformity detected in the phases and phase differences for / = 20 is not particularly significant, 
however, the kurtosis for this distribution is not especially large. 

Phases were drawn from wrapped Cauchy distributions in order to test the method against skies with larger 
values of kurtosis. The wrapped Cauchy distribution has a probability density function given by 



K = -TT L ^- (37) 



© 0000 RAS, MNRAS 000, 000-000 



Phase Correlations in Cosmic Microwave Background Temperature Maps 11 





i 




ii 








iii 








K 


-0.003 


-0.07 


-0.25 


0.06 


0.37 


1.44 


5.76 


oo 


l 
2 


67 


81 


63 


90 


52 


46 


37 


98 


35 


4 


97 


68 


41 


12 


79 


89 


28 


69 


100 


6 


100 


18 


66 


87 


25 


4 


98 


95 


100 


8 


100 


98 


69 


8 


74 


87 


99 


100 


100 


10 


100 


34 


66 


12 


44 


32 


65 


99 


100 


12 


100 


37 


80 


71 


58 


7 


55 


72 


100 


14 


100 


22 


29 


85 


79 


40 


28 


65 


100 


16 


100 


48 


78 


88 


43 


78 


89 


99 


100 


18 


100 


8 


5 


23 


84 


19 


100 


100 


100 


20 


100 


85 


55 


99 


53 


27 


32 


99 


100 



Table 2. Results for non-Gaussian skies based on phase differences between consecutive values of m at a fixed I. The results 
show the significance levels of detected departures from random phases for simulated distributions corresponding to (i) coupled 
phases, (ii) a cardioid distribution and (iii) a wrapped Cauchy distribution. In the latter two cases various parameter choices 
are shown as described in the text. These results arc obtained from 1000 Monte Carlo skies. 



f(8) = tt~ -I. 2 f m ' 0<e<2^,0<p<l (38) 
2ty 1 + p 2 — 2pcos(6) 

The kurtosis is given by 

_ (p 2 - P 4 ) , , 

(i - P y ■ [M) 

As p — > the distribution tends towards a uniform distribution. As p — > 1 the distribution tends towards a distribution 
concentrated at 4> = 0. Fisher (1993) gives a method for simulating the distribution. Sets of coefficients were created 
for p = 0,0.2,0.4,0.6,0.8, 1. The results are given in panel (iii) of Tables 1 and 2. Again p — is not displayed, 
but was used to check the distribution simulation was correct because it should correspond to a uniformly random 
distribution. For p > 0.6 the coefficients are distinctly non-Gaussian. The results in the Table indicate that the 
method can pick out departures from uniformly random phase distributions, at 95% confidence level for \K\ <; 1. 
Note also that the results from p = 1 demonstrate that firm conclusions about the Gaussianity and/or stationarity of 
the sky are difficult to draw using the quadrupole alone. Even for the highly non-Gaussian distribution with K = oo, 
our statistic does not prove to be discriminatory. In summary, however, the tests do show that the method can pick 
out departures from the random phase hypothesis in CMB sky maps using low multipole modes. 



4.2 Quadratic Non-Gaussian Maps 

In order to further examine the performance of our statistic, we applied it to non-Gaussian models with stronger 
physical foundations. At large scales, the CMB is dominated by the Sachs- Wolfe (Sachs & Wolfe 1967) effect so that 
the temperature fluctuations are proportional to fluctuations in the gravitational potential on the last scattering 
surface. Non-Gaussian CMB skies can be generated by introducing a non-linear term in the gravitational potential 
$(n) 

$(n) = *z» + /„,(*£(«)" < *£(n) >) (40) 

where $i(n) refers to the linear part of the potential, f n i is the non- linear parameter controlling the amount of 
non-Gaussianity and the brackets indicate a volume average (Liguori, Matarrese & Moscardini 2003) . We applied our 
method to these types of skies with varying levels of non-Gaussianity controlled by the parameter / ni . 

The results for the most extreme form of non-Gaussianity (/„; = oo) are shown in panel (i) in Tables 3 and 4. 
Note that our statistic is clearly unsuitable for detecting this form of non-Gaussianity. Even for this level of non- 
Gaussianity the method is incapable of finding any deviations of the phases or phase differences from the random 
phase hypothesis. The reason for this is the very simple form of phase association we are testing here does not pick 
up very sensitively the quadratic correlations manifested by the model given in equation (40); see Watts & Coles 
(2003). We discussed this in Section 2.1, in the context of Fourier phases. The reason for this failure is that the 
method we are using is much more sensitive to departures from statistical homogeneity over the celestial sphere than 
it is to statistically homogeneous non-Gaussianity. To put this another way, our method will only be sensitive to 
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Table 3. Results based on phases and phase differences between consecutive values of m at a fixed I. The results show the 
significance levels of detected departures from random phases for (i) the quadratic non-Gaussian map (/„; = oo), (ii) the ILC 
map, (iii) the cleaned TOH map and (iv) the Wiener filtered TOH map. 
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Table 4. Results based on phase differences between consecutive values of I at a fixed m. The results show the significance 
levels of detected departures from random phases for (i) the quadratic non-Gaussian map (/„; = oo), (ii) the ILC map, (iii) the 
cleaned TOH map and (iv) the Wiener filtered TOH map. 



departures from Gaussian behaviour, such as localized edges or blobs, if these features occupy a substantial fraction 
of the celestial sphere. This result will be useful later on, when we interpret the results we obtain from real data. 

It is also worth mentioning that it would not surprising if a large significance level were attached to one of the 
Di(m) values even if the null hypothesis were correct. Since we are looking at 18 values of m when scanning across 
fixed values of m, statistically speaking, one value of m should lead to a probability value larger than 95 per cent. We 
will not try to assess the significance level of any departure from the null hypothesis in the context of this result as 
our motivation for introducing this form of non-Gaussianity is basically to show that our method struggles to detect 
it. Whatever method one uses to combine probabilities across modes would clearly produce a marginal result in this 
case for the probability that there is any departure at all. We will explore this question in more detail in future work 
(Dineen et al., in prep.). 
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Table 5. Results for COBE DMR data based on phases and phase differences between consecutive values of m at a fixed I. 
The method was applied to 5 realisations of the DMR sky in order to demonstrate the stability of the results. 

4.3 The COBE DMR data 

We now discuss the application of our method to real data in order to see if any signs of non-random phases could 
be found in published sky maps. In this subsection we look at the ai, m obtained from COBE-DMR data and in the 
next subsection we apply the method to a; >m obtained from three all-sky maps derived from WMAP data. 

The DMR (Differential Microwave Radiometer) instrument on board the COBE satellite ushered in the era 
of precision cosmology by measuring primordial temperature fluctuations of order AT/T — 1CP 5 in the cosmic 
microwave background (CMB) radiation. The DMR instrument comprised six differential microwave radiometers: 
two nearly independent channels, labelled A and B, at frequencies 31.5, 53 and 90 GHz. The data from all three 
frequencies have been used to calculate the a; jm corresponding to the CMB signal from regions outside the Galactic 
cut. To keep the presentation brief we discuss the behaviour of a; jm for only the even modes. 

The results of the application of the method on the a ijm are displayed in Table 5. The method was applied to five 
realisations (rotations) of the DMR sky in order to examine whether the results were stable and consequently whether 
the conclusions were robust. In the Table it is clear that the large values of \ 2 (and consequently large significance 
levels) do indeed appear stable but lower values (when the phases appear random) the probabilities depend stronly 
upon the realisation. However, we are only interested in significance levels that are larger than 95% and these appear 
stable for all realisations. 

The statistic indicates that both the phases and phase differences of I = 10 mode appear inconsistent with the 
random phase hypothesis. For all 5 realisations, the \ 2 values are larger than 95% of the MC skies. Figure 3 shows the 
distribution of V obtained for realisation 1 from Table 5 compared with the average distribution of the MC skies for 
I — 10 and 12. The distribution of D m (10) clearly differs from that of the random phase MC skies whereas D m (12) 
appears by eye to be drawn form the same distribution. The results in Table 5 also show that the phases obtained 
from the I = 18 also appear non-random, but the confidence of this result is dependent on the realisation. 

The issue of non-Gaussianity in the COBE-DMR data has been controversial (Ferreira, Magueijo & Gorski 1998; 
Bromley & Tegmark 1999; Magueijo 2000), but is likely to involve systematic problems identified in the data set we 
use. Note however that the strong detection we have found is at a different harmonic mode I = 10 than that originally 
claimed by Ferreira et al. (1998) who used a diagnostic related to the bispectrum. There need not be a conflict here 
because, as we discussed above, we are not sensitive to the same form of phase correlations in this method as the 
bispectrum. 

We do not wish to go further into the statistical properties of the COBE-DMR data here. The results we have 
presented are from the original data which is now known to have suffered from some systematic problems (Banday, 
Zaroubi & Gorski 2000). The effect we are seeing is therefore probably not primordial, making this largely an academic 
discussion. The data have also partly been superseded by WMAP. It does, however, illustrate the importance of using 
complementary approaches to identify all possible forms of behaviour in data sets even when they are of experimental, 
rather than cosmological, origin. 

4.4 WMAP results 

The WMAP instrument comprises 10 differencing assemblies (consisting of two radiometers each) measuring over 
5 frequencies (~23, 33, 41, 61 and 94 GHz). The WMAP team have released an internal linear combination (ILC) 
map that combined the five frequency band maps in such a way to maintain unit response to the signal CMB while 
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Figure 3. Examples of the distribution of Kuiper's statistic V for D m (lO) (left) and D m (12) (right) in the COBE-DMR data. 
The red squares and blue stars correspond to the average distribution of V for the 1000 MC skies and the distribution for the 
DMR data respectively. 



minimising the foreground contamination. The construction of this map is described in detail in Bennett et al. (2003). 
To further improve the result, the inner Galactic plane is divided into 11 separate regions and weights determined 
separately. This takes account of the spatial variations in the foreground properties. The final map covers the full-sky 
and should represent only the CMB signal, although it is not anticipated that foreground subtraction is perfect. 

The WMAP ILC map is issued with estimated errors in the map-making technique. This uncertainty could be 
included in the Monte Carlo simulations of the null hypothesis used to construct sampling distributions of our test 
statistic. This would be difficult to do rigorously as errors are highly correlated, but is possible in principle. We have 
decided not to include the experimental errors in this analysis because it is just meant to illustrate the method: a 
more exhaustive search for CMB phase correlations (for which such considerations would be relevant) will have to 
wait for cleaner maps than are available now. 

Following the release of the WMAP 1 yr data Tegmark, de Oliveira-Costa & Hamilton (2003; TOH) have produced 
a cleaned CMB map. They argued that their version contained less contamination outside the Galactic plane compared 
with the internal linear combination map produced by the WMAP team. The five band maps are combined with 
weights depending both on angular scale and on the distance from the Galactic plane. The angular scale dependence 
allows for the way foregrounds are important on large scales whereas detector noise becomes important on smaller 
scales. TOH also produced a Wiener filtered map of the CMB that minimisizes rms errors in the CMB. Features 
with a high signal-to-noise are left unaffected, whereas statistically less significant parts are suppressed. While their 
cleaned map contains residual Galactic fluctations on very small angular scales probed only by the W band, these 
fluctuations vanish in the filtered map. 

The three all-sky maps are available in HEALPix* format (Gorski, Hivon & Wandelt 1999). We derived the ai >m 
for each map using the anafast routine in the HEALPix package. The results are shown in panels (ii) to (iv) in 
Tables 3 and 4. The results for the ILC map suggest there are departures from the random phase hypothesis for the 
phases of / = 12 and 16 and the phase differences for 1=14 and 16. The departure for 1=16 appears very strong with 
the x 2 value being larger than 995 of the values obtained for the MC skies. For these data we have also explored the 
effect of measuring differences between the phases at different I for the same value of m. (We did not present results 
in this case for the COBE-DMR data, as we only discussed the even i-modes in that case and differences between 
adjacent I would require every mode.) Notice that the confidence levels of non-uniformity are larger for D m at fixed 
I than for Di at fixed m. This is an interesting feature of the WMAP data. 

To investigate this mode further we therefore reconstructed the temperature pattern on the sky solely using only 
the spherical harmonic modes for I = 16 with their appropriate ai tm and compared this with a map with the same 
power spectrum but with random phases (figure 4). The result is a striking confirmation of the above argument that 
our method is very sensitive to departures from statistical homogeneity. The / = 16 modes are clearly forming a 
more structured pattern than the random-phase counterpart. The alignment of these modes with the Galactic plane 
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is striking; it may be that the angular scale corresponding to I — 16 may relate to the width of the zones used to 
model the Galactic plane. 

The two maps produced by Tegmark et al. (2003; TOH) also suggest that the phases for the I — 16 scale are 
non-random. The cleaned map shows departures from the random phase hypothesis at greater than 95% confidence 
for the phases of 1=16 and phase differences of I = 14 and 16. The Wiener filtered map finds departures for the 
phases corresponding to I = 16 and phase differences in 1=14. Both the cleaned and the Wiener filtered maps of TOH 
display similar band patterns to that observed in the I = 16 map for the ILC; we show the Wiener filtered TOH map 
for comparison. 

The morphological appearance of these fluctuations as a band across the Galactic plane may be indicative of 
residual foreground still present in the CMB map after subtraction. This is the most likely explanation for the 
alignment of features relative to the Galactic coordinate system, since the scanning noise is not aligned in this way. 
Further evidence for this is provided by Naselsky et al. 2003 from cross-correlations of the phases of CMB maps with 
various foreground components. In any case one should certainly caution against jumping to the conclusion that this 
represents anything primordial. These are preliminary data with complex variations in signal-to-noise across the sky. 
Nevertheless, the results we present here do show that the method we have presented is useful for testing realistic 
observational data for departures from the standard cosmological statistics. 



5 DISCUSSION AND CONCLUSIONS 

In this paper we have presented a new method of testing CMB data for departures from the homogeneous statistical 
behaviour associated with Gaussian random fields generated by inflation. Our method uses some of the information 
contained within the phases of the spherical harmonic modes of the temperature pattern and is designed to be most 
sensitive to departures from stationarity on the celestial sphere. The method is relatively simple to implement, and 
has the additional advantage of being performed entirely in "data space". Confidence levels for departures from 
the null hypothesis are calculated forwards using Monte Carlo simulations rather than by inverting large covariance 
matrices. The method is non-parametric and, as such, makes no particular assumptions about data. The main 
statistical component of our technique is a construction known as Kuiper's statistic, which is a kind of analogue to 
the Kolmogorov-Smirnov test, but for circular variates. 

In order to illustrate the strengths and weaknesses of our approach we have applied it to several test data sets. 
To keep computational costs down, but also to demonstrate the usefulness of non-parametric approaches for small 
data sets, we have concentrated on low-order spherical harmonics. 

We first checked the ability of our method to diagnose non-uniform phases of the type explored by Watts, 
Coles & Melott (2003) and Matsubara (2003). The method is successful even for the low-order modes that have 
few independent phases available. We then applied the method to quadratic non-Gaussianity fields of the form 
discussed by Watts & Coles (2003) and Liguori, Matarrese & Moscardini (2003). Our method is not sensitive to 
the particular form of phase correlation displayed by such fields as, although non-Gaussian, they are statistically 
stationary. Phase correlations are present in such fields, but do not manifest themselves in the simple phase or 
phase-difference distributions discussed in this paper. These two examples demonstrate that our method is a useful 
diagnostic of statistical non-uniformity and its applicability to non-Gaussian fields is likely to be restricted to cases 
where the pattern contains strongly localized features, such as cosmic strings or textures (e.g. Contaldi et al. 2000). 

We next turned our attention to the COBE-DMR data. Claims and counter-claims have already been made 
about the possible non-Gaussian nature of these data, the most likely explanation of the observed features being 
some form of systematic error. Our test does show up non-uniformity in the phase distribution at high confidence 
levels for the I = 10 mode and, less robustly, at / = 18. This is a different signal to that claimed previously to be 
indicative of non-Gaussianity. Its interpretation remains unclear. 

Finally we applied our method to various representations of the WMAP preliminary data release. In this case 
we find clear indications at / = 16 and, less significantly, at I = 14. The distribution of the I = 16 harmonics on the 
sky shows that this result is not a statistical fluke. The data are clearly different from a random-phase realisation 
with the same amplitudes and there is certainly the appearance of some correlation of this pattern with the Galactic 
coordinate system. We are not claiming that this proves there is some form of primordial non-Gaussianity in this 
dataset. Indeed the WMAP data come with a clear warning its noise properties are complex so it would be surprising 
if there were no indications of this in the preliminary data. It will be useful to apply our method to future releases 
of the data in order to see if the non-uniformity of the phase distribution persists as the signal-to-noise improves. 
At the moment, however, the most plausible interpretation of our result is that it represents some kind of Galactic 
contamination, consistent with other (independent) claims (Dineen & Coles 2003; Chiang et al. 2003; Eriksen et al. 
2003; Naselsky et al. 2003). 
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Figure 4. Top: Internal linear combination map constructed with a,[ m for / = 16 only. Middle: The TOH Wiener-filtered map 
with a ; m for I = 16 only. Bottom: A random-phase realization of the map including modes with I = 16 only. 
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As a final comment we mention that the aptitude of our method for detecting spatially localised features (or 
departures from statistical homogeneity generally) suggests that it is may be useful as a diagnostic of the repeating 
fluctuation pattern produced on the CMB sky in cosmological models with compact topologies (Levin, Scannapieco 
& Silk 1998; Scannapieco, Levin & Silk 1999; Rocha et al. 2002); see Levin (2002) for a review. We shall return to 
this issue in future work. 
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